Natural and unnatural parity states of small trapped equal-mass two-component 
Fermi gases at unitarity and fourth-order virial coefficient 
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Equal-mass two-component Fermi gases under spherically symmetric external harmonic confine- 
ment with large s-wave scattering length are considered. Using the stochastic variational approach, 
we determine the lowest 286 and 164 relative eigenenergies of the (2, 2) and (3, 1) systems at uni- 
tarity as a function of the range ro of the underlying two-body potential and extrapolate to the 
ro ^ limit. Our calculations include all states with vanishing and finite angular momentum 
L (and natural and unnatural parity H) with relative energy up to lO.Sftfi, where Q denotes the 
angular trapping frequency of the external confinement. Our extrapolated zero-range energies are 
estimated to have uncertainties of 0.1% or smaller. The (2,2) and (3,1) energies are used to de- 
termine the fourth-order virial coefficient of the trapped unitary two-component Fermi gas in the 
low-temperature regime. Our results are compared with recent predictions for the fourth-order virial 
coefficient of the homogeneous system. We also calculate small portions of the energy spectra of the 
(3, 2) and (4, 1) systems at unitarity. 

PACS numbers: 



I. INTRODUCTION 

Small trapped Fermi gases with contact or short-range 
interactions have attracted a great deal of attention re- 
cently Using lithium or potassium, for example, 
equal-mass two-component systems can be realized ex- 
perimentally by occupying two different hyperfine states. 
For typical experimental conditions, p-wave or higher 
partial wave interactions between two like atoms (say, 
two spin-up atoms) and between two unlike atoms (a 
spin-up and a spin-down atom) are negligibly small. 
Furthermore, by tuning an external magnetic field in 
the vicinity of a Fano-Feshbach resonance, the s-wave 
scattering length Os can be adjusted to essentially any 
value 0. In this paper, we consider the regime where the 
s-wave scattering length is much larger than the range 
To of the underlying two-body model potential. In the 
limit that tq goes to zero and as goes to infinity, the uni- 
tary regime is realized. In this regime, the only mean- 
ingful length scale of the system is given by the oscilla- 
tor length Ohn that characterizes the external confining 
potential [H y|. Throughout, we assume a spherically 
symmetric harmonic potential with angular trapping fre- 
quency f2 (i.e.. Oho = ^Jh/{mn) with m denoting the 
atom mass). 

From a theoretical point of view, small harmonically 
trapped Fermi gases with central short-range interactions 
are particularly appealing since they can be treated with 
comparatively high accuracy by a variety of methods, 
including techniques that have been developed in the 
context of atomic physics, nuclear physics and quantum 
chemistry problems |4l4l6| . For the harmonically trapped 
equal-mass system, the center of mass degrees of free- 
dom separate. Furthermore, the relative orbital angular 
momentum quantum number L, the projection quantum 
number M and the parity 11 are good quantum numbers. 
This implies that the Hilbert space can be divided into 



subspaces, which significantly reduces the complexity of 
the calculations compared to, for example, systems con- 
fined to move within a box with periodic boundary condi- 
tions ^17i |. The harmonically trapped Fermi gas consist- 
ing of two spin-up and two spin-down atoms with vanish- 
ing angular momentum has been treate d by a variety of 
techniques in the literature (see Refs. [UllSlj for reviews). 
The ground state energy and ground state properties of 
the (2, 2) system in the zero-range limit, for example, are 
by now well characterized [E, Much less, however, 
is known about the excitation spectrum [1, [l9l |. 
which contains both natural and unnatural parity states, 
i.e., states with parity 11 = (—1)^ and 11 = (—1)^+^, re- 
spectively. While a good portion of the excitation spec- 
trum of the (2, 2) system with natural parity has been de- 
termined throughout the crossover and at unitarity [T^ . 
little is known about the properties of states with unnat- 
ural parity. Moreover, the energy spectra of the (3,1), 
(3, 2) and (4, 1) systems have not yet been characterized 
in detail. 

This paper presents extensive benchmark results for 
the (2,2) and (3,1) energies of natural and unnatural 
parity states at unitary. In addition, we present portions 
of the energy spectra of the (3, 2) and (4, 1) systems. We 
then use the energy spectra of the (2,2) and (3, 1) sys- 
tems at unitarity to determine the fourth-order virial co- 
efficient 64 of the trapped system in the low-temperature 
regime. The fourth-order virial coefficient enters into the 
virial equation of state, which allows for the determina- 
tion of the tmiversal thermodynamics of two-component 
Fermi gases in the temperature regime down to about 
half the Fermi temperature Tp |20i - (2^ . For the tem- 
perature regime in which we have convergence, i.e., for 
ksT < 2Ml/i, where T denotes the temperature and 
kB Boltzmann's constant, we find that the fourth-order 
virial coefficient 64 of the trapped system is negative 
and decreases monotonically with increasing tempera- 
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ture. If we assume that 64 continues to change mono- 
tonically with increasing temperature in the medium- 
and high-temperature regime, our resuhs predict that 
the fourth-order virial coefficient of the trapped system 
and — through apphcation of the local density approxima- 
tion (LDA) — that of the homogeneous system approach 
a negative value in the high-temperature limit. This is 
in contrast to recent results (26-28] based on the equa- 
tion of state, determined both experimentally and calcu- 
lated via a diagrammatic Monte Carlo technique. These 
studies predict that the fourth-order virial coefficient of 
the homogeneous system is positive. The discrepancy 
would be resolved if the fourth-order virial coefficient 
of the trapped system was changing non-monotonically 
with temperature, allowing for a sign change of 64 in 
the medium- or high-temperature regime. Analogous 
non-monotonic behavior was found for one of the third- 
order virial coefficients of the trapped unequal-mass two- 
component Fermi gas at unitarity [22]. While we do not 
have access to sufficiently large portions of the energy 
spectra of the (2, 2) and (3, 1) systems to determine the 
fourth-order virial coefficient of the trapped system in 
the medium- and high-temperature regimes (thereby pre- 
venting us from drawing definite conclusions) , our results 
illuminate a number of aspects related to the determina- 
tion of the virial coefficients from few-body energy spec- 
tra. 

Section In] introduces the system under study, reviews 
the stochastic variational approach, and presents details 
regarding our implementation. Compact expressions for 
the relevant matrix elements for natural and unnatural 
parity states are presented in Appendix [Xj Section IIIII 
(see also supplementary material [30]) summarizes our 
extrapolated zero-range energies for the (2,2), (3,1), 
(3,2) and (4,1) systems in tabular form and discusses 
their characteristics. Section HVl uses the (2, 2) and (3, 1) 
energies to determine the fourth-order virial coefficient of 
the trapped system at unitarity. Lastly, Sec. IVl concludes. 



II. SYSTEM UNDER STUDY AND 
STOCHASTIC VARIATIONAL APPROACH 

We consider a two-component Fermi gas with ni spin- 
up and 712 spin-down atoms of mass m with n = ni+n2. 
We assume that the atoms are confined by a spherically 
symmetric trapping potential with angular frequency fl. 
Furthermore, we assume that the spin-up and spin-down 
atoms interact through a short-range interaction poten- 
tial Vth{rpq), where rp (p = 1, ■ ■ ■ ,n) denotes the position 
vector of the p*'' atom measured relative to the center of 



where 



the trap and 



Tq], and that atoms with like 



spins do not interact. The model Hamiltonian H then 
reads 



H 



E 

p=i 



1 



2m 



(1) 



(2) 



p=l q=ni+l 



Throughout, we are interested in the regime where the 
s-wave scattering length as of the interspecies inter- 
action potential Vtb becomes infinitely large. For the 
(rii,ri2) — (1,1) and (2,1) systems, semi-analytical so- 
lutions are known if Vth coincides with the zero-range 
(5-function potential [1, |3l|. For (ni,n2) systems with 
ni + 712 > 4, however, no such semi- analytical solutions 
are known. To determine the eigenenergies of (741,712) 
systems with ni+n2 = 4 and 5, we separate off the center 
of mass motion and resort to a numerical technique, the 
stochastic variational approach [32| . In this approach, it 
is convenient to model the interactions between the un- 
like atoms through a Gaussian potential Vg (r) with depth 
-Vq {Vq > 0) and range rg [1], 



^gi'^) = -Vbexp 



V2; 



ro 



(3) 



To treat the unitary system, we adjust the depth Vq of 
Vg for a given rp such that the two-body system in free 
space supports one zero-energy s-wave bound state but 
no deep-lying bound states. To determine the zero-range 
energies, we consider a number of rg, rg ^ Oho, and 
extrapolate the finite-range energies to the rg — >■ limit 
(see Sec. IIIII for examples). 

We take advantage of the fact that the Hamiltonian H 
separates into the center of mass Hamiltonian H"'^^ and 
the relative Hamiltonian H"'^ H = H"'^ + In the 

following, we consider the relative Hamiltonian and 
use the stochastic variational approach to determine the 
eigenenergies and eigenstates of the Schrodinger equation 
^roi^rci ^ ^rci^^_^^roi^ ^{eTe, wc exphcitly indicate the 
dependence of the eigenenergies on ni and 712 but, for 
notational simplicity, not that of the Hamiltonian and 
the wave function. To compact the notation, we write 
H'"^ as H'"^ = T'"^ + V^,% + Fint, where T"""" denotes the 
kinetic energy operator associated with the relative mo- 
tion, and Vj^ap the contribution of the confining potential 
associated with the relative degrees of freedom. 

The stochastic variational approach is a basis set ex- 
pansion approach that writes the relative wave function 
of a given state in terms of a set of basis functions 



rcl 



E 

fe=i 



(4) 



Here, the Cfc denote expansion coefficients and A an anti- 
symmetrization operator that ensures that the wave func- 
tion is anti-symmetric under the exchange of any pair of 
like fermions. In Eq. ([3]), N}, denotes the number of basis 
functions. As with other basis set expansion approaches. 
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the Ritz variational principle ensures that the lowest en- 
ergy as well as the higher-lying energies obtained by 
the stochastic variational approach are rigorous upper 
bounds to the exact eigenenergies of the system [3^ . In 
the following, we introduce the basis functions used in 
this work, which have good orbital angular momentum 
L, projection quantum number M and parity 11; here, L, 
M and 11 are associated with the relative motion. 

Following Refs. [33 - l38| . we write the basis functions 
i>k as a product of a correlated Gaussian [second line 
of Eq. ©] and a "prefactor" [first line of Eq. (O] that 
carries the angular momentum L of the system, 

Mx) = \vik\^'\v2kt'[Yi,{vik)®Yi,{v2k)]LM 

xexp(-:^). (5) 

Here, x collectively denotes the n — 1 Jacobi vectors 
Pp, where p = I,-- - ,n — 1. The notation \Yi^{vik) ® 
Yi2{^2k)]LM indicates that the spherical harmonics Yi^mi 
and Yi^m2 ^re coupled to form a function with angular 
momentum L and projection quantum number M. For 
states with natural parity, i.e., for states whose parity is 
given by n = (—1)^, we choose li = L and I2 — Q |32| - 
|35|. For states with unnatural parity [L > 0), i.e., for 
states whose parity is given by 11 = (—1)^+^, we choose 
li = L and ^2 = 1 [3^.l37t. The basis functions that de- 
scribe unnatural parity states with L = have a slightly 
different form since the construction of states with L — Q 
and n = — 1 requires the coupling of three spherical har- 
monics with li, I2 and ^3 > [Sa,!!!]. The matrix Af. 
is symmetric and positive-definite, and has dimensions 
(n — 1) X (n — 1). The n(n — l)/2 independent elements 
of Af^ are treated as variational parameters and opti- 
mized semi-stochastically. The three-dimensional vectors 
vik and V2k, referred to as global vectors since they de- 
pend on all n — 1 Jacobi vectors, are defined through 
vik = Z]p=i "ife,pPp = ^ifc^ aiid similarly for V2k- The 
vectors uik and U2k are optimized semi-stochastically, 
where uik = (wifc,i, • • • , uife,„-i) and similarly for U2k- 

A key benefit of the basis functions given in Eq. ([S]) is 
that the overlap matrix element Ok'k = {'i'k' iV'fe), the ma- 
trix element for the kinetic energy operator (T''°')fe'fc = 
('i/'fe' |T'''''|'!/;fe), and the matrix element for the confining 
potential (K^!,p)fc'fc — ("^fc' l^traplV'fc) reduce to compact 
expressions I32}437|. Here, it is understood that the in- 
tegration is performed over all 3n — 3 Jacobi coordinates 
and that i/'fe is characterized by A^., uik and U2k while if^k' 
is characterized by Ay, uik' and U2k'- Moreover, a com- 
pact expression can also be found for the matrix elements 
(Mnt)fc'fc = {ij-'k' iVintlipk) associatcd with the atom-atom 
interaction if Vth is modeled by the Gaussian potential 
Vg . Appendix |X] summarizes explicit expressions of the 
matrix elements with natural parity (any L) and unnat- 
ural parity (L > 0). The matrix elements for states with 
0~ symmetry can be found in Refs. [36l. Issj. 

We note that the overlap matrix element Ok'k between 
two different basis functions does not vanish, i.e., the ba- 
sis set employed is not orthogonal. This implies that the 



determination of the eigenenergies amounts to the diag- 
onalization of a generalized eigenvalue problem defined 
by the Hamiltonian and overlap matrices [s^l • While one 
might think, at first sight, that the non-orthogonality of 
the basis functions could introduce numerical instabili- 
ties, it has been shown in previous work [^, [^, [l^, 
that numerical instabilities due to linear dependence is- 
sues can be avoided completely for the systems of interest 
in this work if the basis sets are chosen carefully. 

Our strategy to optimize the large number of non- 
linear variational parameters is quite simple 32]. We 
start with a reference basis set, which could consist of 
just one basis function or as many as several 100 or 1000 
basis functions. We then enlarge this reference basis set 
by one basis function, which is chosen from a large num- 
ber of trial basis functions, typically between several 100 
and several 1000. Each trial function is characterized by 
a different set of variational parameters. To decide which 
trial basis function to keep, we calculate the energy for 
each of the enlarged trial basis sets, which consist of the 
reference basis set plus one of the trial basis functions, 
and choose the one that results in the largest reduction 
of the energy of the state of interest. The state of interest 
could be the ground state or an excited state. The proce- 
dure is repeated till the basis set is sufficiently complete 
to describe the state of interest with the desired accuracy. 

When optimizing a state whose energy is nearly de- 
generate with that of another state or when optimizing 
highly excited states, some care needs to be exercised. 
In the former case, we find it advantageous to optimize 
two or more states simultaneously. In the latter case, we 
find it beneficial to start with a basis set that provides 
a reasonably accurate description of the lower lying part 
of the energy spectrum. The advantage of our optimiza- 
tion procedure is that the basis set is optimized for a 
particular state or a particular subset of states. Corre- 
spondingly, we work with comparatively small basis sets. 
The energies of the (2,2) and (3,1) systems at unitar- 
ity (see Table |T] and supplemental material) are obtained 
using basis sets that consist of 700-3400 basis functions, 
while the energies of the (3, 2) and (4, 1) systems (see 
supplemental material) are obtained using basis sets that 
consist of 1500-3800 basis functions. 



III. ENERGIES OF SMALL TRAPPED FERMI 
GASES 

One key purpose of this paper is to elucidate how 
we determine a large portion of the energy spectrum of 
trapped two-component Fermi systems with n = 4 and 
5, and to tabulate the extrapolated zero-range energies. 
We believe that the tabulation of the energies is useful 
as these energies provide much needed highly accurate 
benchmark results that can be used to assess the accu- 
racy and validity regime of alternative approaches. We 
anticipate that the tabulated energies will also prove use- 
ful in other applications. 
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FIG. 1: (Color online) Illustration of convergence for the (3, 1) 
system at unitarity with 1"*" symmetry and ro = 0.04aho. 
Solid and dashed lines show the quantity Ae'^^ , where Ae^^^i = 
[E'a^iiNb) - El';i{Nb oo)]/E'3'^{{Nb ->■ oo), for states 1 and 
12 as a function of l/Nt- Dotted lines show the extrapolation 
to the A^b — !■ oo limit. The inset shows a blow-up of the small 
1/Nb region. 



Figure [T] shows an example of our basis set opti- 
mization for the (3,1) system with 1"*" symmetry and 
ro = 0.04aho at unitarity. SoUd and dashed hnes show 
the fractional difference Ae™l for the ground state (state 
1) and state 12 (s^l, respectively, between the relative 
energy E'™} for a basis set of size iV^ and the energy for 
an infinite basis set. The dotted lines in Fig. [T] show the 
extrapolation to the A^^ — )■ oo limit. It can be seen that 
the ground state energy converges notably faster than 
the excited state energy. The energies for iVf, — 800 and 
Nb = 900 are El'^\{ro = 0.04aho) = 5.08294:hn for state 
1 and E'^'^\{ro = 0.04aho) = 10.1788?ir2 for state 12, re- 
spectively. The basis set errors for these basis sizes are 
0.0002% and 0.003%, respectively, i.e., the energies of 
states 1 and 12 he respectively O.OOOOl^ifi and 0.0003;if7 
above the extrapolated energies for the infinite basis set. 
The low-lying states of the (3,1) system with 1+ symme- 
try at unitarity converge relatively quickly with increas- 
ing Nf,. The convergence is slower for most other states 
and, in general, we choose the size of our basis sets for 
the (2, 2) and (3, 1) systems such that the basis set ex- 
trapolation error is smaller than 0.1 %. 

Figure [5] exemplarily illustrates the range dependence 
for the relative energy of the (3, 1) system with 1+ 
symmetry. Figure HJa) shows the range dependence of 
the ground state energy. Fig. ^h) shows the range de- 
pendence of the energy associated with state 12, and 
Fig. [He) shows the range dependence of the energy 
for a state that depends comparatively weakly on rg 
(state 5). In Figs. OJa) and^b), the energies vary to a 
very good approximation linearly with vq for sufficiently 
small fp /o-hn- This finding is in agreement with earlier 
work [alllMlIl- For the ground state [see Fig.dja)], 
the range dependence is quite weak and linear behavior 
is only observed for tq < O.OSoho- 
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FIG. 2: (Color online) Illustration of finite-range dependence 
for the (3, 1) system with 1"*" symmetry at unitarity. Squares 
show the relative eigenenergies E^'^{{Ni,) for various ranges ro 
of the underlying two-body interaction potential for (a) the 
ground state (state 1), (b) state 12, and (c) state 5; A^;, is the 
largest basis set considered. The energies provide variational 
upper bounds and the estimated basis set extrapolation error 
is indicated by errorbars; in panels (a) and (b), the basis set 
extrapolation error is smaller than the symbol size and thus 
not visible. In panels (a) and (b), solid lines show linear fits 
to the energies El'^\{Nb) [the fit shown in panel (a) includes 
the energies for the five smallest ro values]. 



In Fig. [Ifc), the zero-range energy agrees to within 
0.00002ftri with the energy of the non-interacting sys- 
tem. This, combined with the very weak dependence of 
the energy on ro and the fact that the energy approaches 
the zero-range limit from below, suggests that this state 
is not affected by s-wave scattering but only by higher- 
partial wave scattering. In the zero-range limit, energy 
shifts associated with higher-partial wave scattering pro- 
cesses vanish. Our interpretation is corrobated by a per- 
turbative calculation along the lines of that performed 
in Refs. [1, [l2| , which utilizes zero-range contact interac- 
tions. For the (3, 1) system with = 1+ symmetry, we 
find, in agreement with our results based on the stochas- 
tic variational approach, that there exists one state with 
relative energy 17hfl/2 and six states with relative energy 
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21KI/2 that are independent of as- 

We refer to states that are unaffected by s-wave inter- 
actions as unshifted states. We find that a relatively large 
number of states fall into this category. Their existence 
and likelihood of occurance has already been discussed 
for the (2, 1) and (2,2) systems in the literature [3. [l^. 
For the (2, 1) system, e.g., all unnatural parity states are 
unaffected by s-wave interactions in the zero-range limit. 
For the (2, 2) and (3, 1) systems, unnatural parity states 
can be affected by s-wave interactions in the zero-range 
limit. The only exception are states with 0~ symmetry, 
which are unshifted. This behavior can be intuitively un- 
derstood within a picture that utilizes angular momen- 
tum coupling. To construct a state with 0~ symmetry, 
the coupling of three finite angular momenta is needed. 
These angular momenta can be envisioned as being each 
associated with one of the three Jacobi vectors that char- 
acterize the n = 4 system. As a consequence, the s-wave 
interactions are effectively turned off by the nodal struc- 
ture of the wave function. For n = 5, this argument 
predicts that states with 0~ symmetry can be affected 
by s-wave interactions since the system is characterized 
by one more Jacobi vector than angular momenta needed 
to ensure the 0~ symmetry. Indeed, this prediction is in 
agreement with our results from the perturbative and 
stochastic variational calculations. 

Table U summarizes our extrapolated zero-range ener- 
gies Ef\{ro = 0), Ef\{rQ = 0) < lO.S^rj, for states with 
1+ symmetry at unitarity that are affected by s-wave in- 
teractions. The zero-range energies are obtained by cal- 
culating the energies of a given state for several ranges tq 
between 0.0025 < ro/flho < 0.08 and by then fitting these 
energies for the largest basis set considered by a linear 
function. The third column in Table U shows the slopes 
X, which characterize the dependence of E^^\ on ro at 
unitarity. We find that the slopes for states that are af- 
fected by s-wave interactions are positive. Table U shows 
that the slopes vary over nearly two orders of magnitude. 
The slopes can be related to the effective range Teff using 
the relation Toff = 2.032ro. This numerically determined 
relationship is specific to the Gaussian model potential 
employed in this paper and is quite accurate over the rg 
values considered. It may be used to estimate the leading 
order dependence of the energies on the effective range 
for the Gaussian model potential. 

The relative energies at unitarity for zero-range inter- 
actions can be written in the form {2q-\-SL,v + ^)fi'^ jl.l40|. 
where slm is associated with the eigenvalue of the hyper- 
angular Schrodinger equation and where the radial quan- 
tum number q takes the values 0, 1, • • • (although the sl,i/ 
depend on 11, this dependence is not explicitly indicated 
for notational simplicity) . The fourth column of Table H] 
shows the sl,v values determined from our energies for 
q = 0, i.e., for the lowest rung of the ladder with 2qMl 
spacings. The extrapolated zero-range energies of states 
2 and 7, e.g., lie 2mQlMl and 4.0006^1], respectively, 
above the energy of the ground state. Correpondingly, 
we assign the quantum numbers q = 1 and g = 2 to these 



TABLE I: Relative energies Ef], for the (3, 1) system with 
= 1"*" symmetry [only states that are affected by s-wave 
interactions are included; each energy is (2L -I- l)-fold de- 
generate]. The first column indicates the state number (st. 
no.). The second column shows the extrapolated zero-range 
energy -E3^i(ro — 0) at unitarity; the uncertainty is estimated 
to be 0.1 % or smaller. The third column indicates the de- 
pendence of the energy at unitarity on the range ro of the 
Gaussian potential Vg. We assume a linear dependence and 
write El%ro) = -E^°i(ro = 0) x(^o/aho)ftn. The fourth col- 
umn shows the sl,v value determined from the energy; the 
value of sl,v is only shown for the lowest rung of a ladder, 
i.e., for states with g = 0. The last column shows sj,'^ of 
the non-interacting state that is "paired" with the interact- 
ing state when determining AQa.i (see Sec. lIVp . There exist 
1 and 6 unshifted states with energy nh£l/2 and 21/if2/2, 
respectively. 
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5.0819 


0.04 


4.0819 


5.5 
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7.0820 


0.03 
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7.6056 


0.51 


6.6056 


7.5 
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8.1456 


0.76 


7.1456 


7.5 
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8.9846 


1.19 


7.9846 


9.5 
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9.0825 


0.03 
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9.1324 


0.28 


8.1324 


9.5 


9 


9.4544 


0.46 


8.4544 


9.5 


10 


9.6060 


0.55 






11 


9.6847 


1.17 


8.6847 


9.5 


12 


10.147 


0.80 







states, i.e., we identify them as belonging to the same 
ladder as the ground state. The small deviations from 
the 2qMl spacings can be interpreted as a measure of our 
numerical accuracy. For the states considered in Table U 
the 2qhVt spacing is fulfilled to better than 0.1 %. We 
find that the energies of states that belong to the same 
ladder are characterized by similar slopes. 

For some symmetries, nearly degenerate states exist 
in the energy range £'''°' < IQ.bKl. Figure |3] shows the 
range dependence of the (3,1) energies with 3^ symmetry 
corresponding to states 15 and 16. This figure illustrates 
exemplarily that the "ordering" of states can change as a 
function of rg, i.e., that the energies of two or more states 
can cross at finite rg. Crossings like these can only be 
resolved by considering at least three different tq values 
for each state. 

Following the format of Table HI the supplemental ma- 
terial tabulates the energies of the (2,2) and (3, 1) sys- 
tems. The results are obtained by analyzing the finite- 
range energies determined by the stochastic variational 
approach along the lines discussed above. For the (2, 2) 
and (3, 1) systems, there exist 286 and 164 states at 
unitarity with relative energy E'^^ smaller or equal 
to lO.Bhn [not counting the (2L + 1) multiplicity]. Of 
these states, respectively 52 and 46 are unshifted. The 
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FIG. 3: (Color online) Illustration of finite-range dependence 
for (3, 1) system with 3~ symmetry at unitarity. Circles and 
squares show the relative eigenenergies El';i{Nb) for states 
15 and 16, respectively, for three different ranges ro of the 
underlying two-body interaction potential; A'^i, is the largest 
basis set considered (typically, A^i, increases with decreasing 
ro/a-ho)- The energies provide variational upper bounds and 
the estimated basis set extrapolation error is indicated by 
errorbars. Solid lines show linear fits to the energies El°\{Ni,). 
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shifted energies are characterized by respectively 170 and 
89 sl,l> values. FiguresSfa) andllfb) show the density of 
states of the (2, 2) system and the (3, 1) system, respec- 
tively, at unitarity. The plots account for the (2L -I- 1)- 
multiplicity, and the density of states is shown separately 
for the shifted and unshifted states. It can be seen that 
the density of states increases significantly with increas- 
ing energy for both the (2, 2) and (3, 1) systems. 

The supplemental material also tabulates results for 
the (3, 2) and (4, 1) systems. For these systems, the 
convergence is slower than for the n = 4 systems, and 
we choose the size of our basis sets such that the basis 
set extrapolation error is smaller than 1 %. Since the 
calculations for n = 5 are significantly more demand- 
ing than for n = 4, we restrict ourselves to states with 
Ktn2 - l7hn/2. We first extrapolate the energy for 
a given tq to the infinite basis set limit, and then de- 
termine the zero-range energy from these extrapolated 
energies. For the (3, 2) and (4, 1) systems, there exist 19 
and 4 states with energies -E^^/n^ ^ 17?if2/2 at unitar- 
ity [not counting the (2L-f 1) multiplicity and excluding, 
for technical reasons, (3,2) states with 0~ symmetry]. 
All of these energies correspond to shifted states. One of 
the (3, 2) energies corresponds to a "repeated state" with 
hyperradial quantum number q = 1. 



IV. 4" -ORDER VIRIAL COEFFICIENT 

This section uses the (2, 2) and (3, 1) energies to deter- 
mine the fourth-order virial coefficient 64 of the s-wave 
interacting two-component Fermi gas under spherically 
symmetric harmonic confinement at unitarity in the low- 



FIG. 4: (Color online) Panels (a) and (b) show the density of 
states for the (2, 2) and (3, 1) systems at unitarity; only the 
relative degrees of freedom are accounted for. The histograms 
show the number of energies corresponding to shifted states 
per hn/4: while the crosses show the number of energies corre- 
sponding to unshifted states. The histograms and the crosses 
account for the (2L -|- l)-multiplicity of the energies. 



temperature regime. We also summarize a few results for 
the low-temperature behavior of the higher-order virial 
coefficients. 

The virial coefficients &„ enter into the virial equation 
of state, which describes the finite temperature behavior 
of trapped two-component Fermi gases |20l- |27l . l29l . l4l1 - i43| . 
We work in the grand canonical ensemble and denote 
the fugacities of component 1 and component 2 by zi 
and Z2, respectively, where the Zi are defined in terms of 
the chemical potentials /i^ of the ith component and the 
temperature T, 

= exp[^i/(fc_Br)]. (6) 

The thermodynamic potential fi^^-' of the harmonically 
trapped Fermi gas can be written in terms of the thermo- 
dynamic potentials fjj^'' and fij^^ of the non-interacting 
components 1 and 2, and an "interaction piece" AQ^^^^ 
that accounts for the interactions between the atoms of 
component 1 and the atoms of component 2 [2^, [2^, l4l| - 

(00 00 \ 
E E bn^.n^^r^r]- (7) 
"1=1 "2 = 1 / 
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TABLE II: The second and third columns show the expres- 
sions for QiAbn and b'j^^ , n = 2 — 5, for the trapped system. In 
deriving these expressions, we used that AQ„j,„2 = A(5„2,„^ 
for the systems considered in this paper. 



n 


QiA6„ 




2 


A0i,i/2 





3 


AQ2,i 


-2fe2Ql 


4 


AQs.i + AQ2,2/2 


-HiQif + b'zQi + 2(52] - 2b3Qi 


5 


A(54,l + AQ3,2 


-262[fe2(Ol)' + Ql02+Q3] 






-63[(Oi)'' + 2Q2 + 262Q1] - 264O1 




In Eq. ([7]), Qi denotes the canonical partition function 
of a single particle in a spherically symmetric harmonic 
trap with angular frequency $7, 

Qi=e3-/2(e--l)-3, (8) 

where w denotes a dimensionless inverse temperature, 

(9) 



If we restrict ourselves to spin-balanced systems with 
equal masses, the fugacities zi and Z2 are equal, z = 
z\ = Z2, and Eq. ([7]) reduces to 



(10) 



where 62 = ^i,i/2, ^3 = (^1,2 + &2,i)/2, ^4 = (^1,3 + &3,i + 
^2,2)72, and so on. 

We find it convenient to write the virial coefficients b„ 



FIG. 5: (Color online) Virial coefficient &2 of the trapped two- 
component Fermi gas at unitarity as a function of the inverse 
temperature a). The solid line shows 62, Eq. (|17|l . while the 
dash-dot-dotted, dash-dotted and dashed lines show 62 ob- 
tained by setting gmax in Eq. (|16|) to 0, 1 and 10, respectively. 
The solid horizontal line shows the high-temperature limit 



are determined by the total energies E!^^^'^^ and -EJJ''^ 
of the interacting two-component and non-interacting 
single-component systems, respectively. It is important 
to note that the energies E!^^'^^ and E^^_^^ contain the cen- 
ter of mass energy. The summation over j in Eqs. (jl3[) 
and extends over all states allowed by symmetry. 
For ni = 1, the sum in Eq. (jl4p can be performed ana- 
lytically, yielding Eq. ([5]). In the high-temperature limit, 
one finds for systems with zero-range interactions at uni- 
tarity that [23] 



'OJ 



(15) 



Abr, 



(11) 



where b^f is determined by the virial coefficients bj and 
the canonical partition functions Qj with j < n. The in- 
teraction piece Abn, in contrast, accounts for the "new" 
physics introduced by the interacting (ni,n2) clusters 
with n = ni + n2. Explicit expressions for A6„ and 
b^^^ are given in Table HIl where the AQn^^n^ are defined 
in terms of the canonical partition functions QJf ' and 
Qm of the interacting (711,712) system and the single- 
component system with ni atoms, respectively. 



AQn,, 



?int 
ni ,712 



(12) 



The temperature-dependent canonical partition func- 
tions Q'^l^^ and Q„i , 



QnU-E^^p[-<"i:4/(^BT)] 

j 



and 



E 

3 



exp[-<'V(A:BT)], 



(13) 



(14) 



The second-order virial coefficient of the trapped sys- 
tem at unitarity takes the simple form [23| 



9max ^ 

b2 = Ab2= Hm V-[ 



q=0 



or, performing the infinite sum. 



b2 = \e-^''{l 



-(2g+l/2)LD 



'(2<j+3/2)ii 



0- 



(16) 



(17) 



The solid line in Fig. [5] shows the second-order virial co- 
efficient 62, Eq. (|17l) . as a function of w. In the high- 
temperature (small w) limit, 62 approaches the constant 

62°'' , b^2^ = 1/4 (solid horizontal line in Fig. [5]), which 
can be obtained by Taylor-expanding Eq. (fT7|) . To illus- 
trate the convergence of 62 with increasing energy cutoff, 
dash-dot-dotted, dash-dotted and dashed lines show 62 
obtained by setting gmax in Eq. (fTB]) to 0, 1 and 10, re- 
spectively. For a finite energy cutoff, it can be seen that 
62 goes to in the small uj region as opposed to 62°'' = 1 /4. 
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FIG. 6: (Color online) Virial coefficient 63 of the trapped 
two-component Fermi gas at unitarity as a function of the in- 
verse temperature iD. The solid line shows 63 with Lmax and 
t'max set to very large values (see Ref. [gsl] for details) while 
the dash-dot-dotted, dash-dotted and dashed lines show bs 
obtained by limiting Lmax and fmax in Eq. 1)18^ such that 
sl,i' < 11/2, < 19/2 and < 50, respectively. The solid hori- 
zontal line shows the high-temperature limit b'^^\ Eq. (|19|l . 



As expected, a larger energy cutoff provides an accurate 
description of 62 over a larger temperature range, i.e., 
down to a smaller inverse temperature Cj. 

The relative three-body energies at unitarity and for 
vanishing s-wave scattering length Os can be written as 
{2q + SL,^ + (see Sec. ^ and {2q + sl\^ + l)hn, 

respectively. Performing the sum over q analytically, the 
interaction piece A&3 of the trapped three-body system 
at unitarity takes the form [2^, [20| 



A6; 



i3 = lim e2"(e2'^ - x 

I-'niax , -t' m ax 00 



(18) 



v=0 L=0 



Using large Lmax and ^'max, a fully converged pointwise 
representation of 63 is obtained [23, (see solid line in 
Fig. in]). Using the analytical forms for Qi and 62, Eqs. (|H]) 
and pri) . we find that 6™^ diverges as —u)~^ /2 + oj~^ /8 in 
the high-temperature limit. This divergence is cancelled 
by a divergence of A63 of opposite sign. As a result, 63 
is well behaved in the small w (high T) limit. A careful 
analysis of the high-temperature behavior gives [2^ 



bf'> = -0.0683396093112849(1) 



(19) 



(see horizontal solid line in Fig. 

To illustrate the convergence of 63 with increasing L^ax 
and fmax [see Eq. (IT^ ]. dash-dot-dotted, dash-dotted 
and dashed lines in Fig. [S] show 63 calculated using A63 
from Eq. (|T8)) with Lmax and fmax chosen such that 
sl,u < 11/5, < 19/5 and < 50, respectively. No cutoff is 
imposed in calculating 6™^. In these calculations, we in- 
clude the S8-iri6 number of sj^ u and ^ in evaluating A63, 
i.e., each, interacting sl v 

value is paired with the corre- 
sponding non-interacting s™ ^ value. Figure [S] shows that 



the cutoff introduces a divergence in 63. This divergence 
arises because the cutoff alters the high-temperature be- 
havior of A63, which implies that the divergencies of 63''^ 
and A63 no longer cancel. Importantly, 63 is converged in 
the low-temperature (large tj) regime even for a relatively 
small cutoff. This allows us to use the converged low- 
temperature tail to constrain 63 in the high-temperature 
regime. Extrapolating 63 (calculated using a cutoff of 9) 

to the high-temperature limit, we find 63°'' ~ —0.068(1), 
which deviates by less than 2% from the exact value. 
The validity of the employed extrapolation scheme cru- 
cially hinges on the fact that the functional form of 63 
changes "predictably" as a) changes from the low- to the 
medium- to the high-temperature regime. For example, 
if 63 changed sign in the medium- or high-temperature 
regime, as is the case for the coefficient 62.1 that char- 
acterizes the behavior of two identical fermions and one 
lighter fermion (with a mass ratio from 3.11 to 8.62) [29| . 
the extrapolation employed above would predict the in- 
correct high-temperature limit of 62,1. 

The interaction piece A64 of the fourth-order virial co- 
efficient can be expressed analogously to A63. In partic- 
ular, we write the energies at unitarity in terms of the 
slm (see Sec. IIIII and the supplemental material for a 
listing of the Sh^v values) and perform, as in the three- 
body case above, the sum over the hyperradial quantum 
number q analytically. Since both natural and unnat- 
ural parity states of the four-body systems are affected 
by the s-wave interactions, the Sh,v values corresponding 
to both natural and unnatural parity states need to be 
included when evaluating A64. The reference piece 64®* 
diverges as 



Lrcf _ l ~ -6 
O4 — —a; 



149 _2 

-bJ 



16 
1 



1 + 64&: 



(0) 



32 



3 c:;-3 



646!,"^ - 5126;, _i 
—a) 



u(2) 



3840 



256 



(20) 



in the high-temperature limit. This divergence must be 
cancelled by a divergence of A64 of opposite sign. 

Dash-dot-dotted, dash-dotted and dashed lines in 
Fig. [7] show 64 at unitarity obtained by using the full 
expression for and by limiting the sums over v and 
L in A64 such that sl,v < 11/2, < 15/2 and < 19/2, re- 
spectively. For the largest cutoff, our calculation includes 
169 and 89 sl,^ values associated with shifted states [not 
counting the (2L -I- l)-multiplicity] of the harmonically 
trapped (2, 2) and (3, 1) systems with zero-range inter- 
actions, respectively. Figure [7] shows that 64 is negative 
in the low-temperature (large uj) regime and that neither 
&4 nor its first or second derivatives with respect to uj 
change sign in the regime where 64 is converged. This 
motivates us to extrapolate the converged part of 64 to 
the medium- and high-temperature regime (see dotted 
line in Fig.[7l), yielding 5^^ = -0.0020(5). 

The LDA predicts that the virial coefficient b^°™ of the 
homogeneous system is related to the high-temperature 
limit of the nth order virial coefficient of the trapped 
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FIG. 7: (Color online) Virial coefficient 64 of the trapped 
two-component Fermi gas at unitarity as a function of the 
inverse temperature ui. The dash-dot-dotted, dash-dotted and 
dashed lines show &4 obtained by limiting sl,v to be smaller 
than 11/2, 15/2 and 19/2, respectively. The dotted line shows 
our attempt to extrapolate to the high-temperature limit; this 
extrapolation assumes that 64 changes "predictably" from the 
low- to the medium- to the high-temperature regime. The 
inset shows the same data as the main figure. In addition, 
the solid horizontal line shows the high-temperature limit b'f^ 
determined by applying the LDA to the fourth-order virial 
coefficient predicted for the homogeneous system [2^ . 



system via [23 



jhom^^3/2W0)^ 



(21) 



Application to our extrapolated 64"'' yields 64°™ 
—0.016(4). This value for the homogeneous system dif- 
fers in both sign and magnitude from the values 64°"^ = 
-1-0.096(15) dl and 6^™ = -h0.065(10) |3 determined 
from experimental data. These experimental values have 
been found to be consistent with the equation of state de- 
termined by a diagrammatic path integral Monte Carlo 
approach [28| . Given the disagreement between our value 
and that reported in the literature, we speculate that 
the fourth-order virial coefficient of the trapped system 
changes sign in the medium- or high-temperature limit, 
implying that the applied extrapolation scheme does not 
predict the correct medium- and/or high-temperature be- 
havior of 64. If this conclusion is correct, it would fol- 
low that the determination of the medium- and high- 
temperature behavior of the fourth-order virial coeffi- 
cient of the trapped systems requires, if determined via 
the microscopic energy spectra, knowledge of large por- 
tions of the energy spectra of the (2, 2) and (3, 1) sys- 
tems. This suggests that other approaches, based on 
Feynman diagrams or based on simulating the finite tem- 
perature behavior directly numerically, may be more suit- 
able than the approach pursued here for determining the 
temperature-dependence of 64. 

We also analyzed the low-temperature tail of 65. Us- 
ing the (3, 2) and (4, 1) energies from the supplemen- 
tal material, we find that the fifth-order virial coeffi- 
cient of the trapped Fermi gas at unitarity is positive 



in the low-temperature limit. High precision measure- 
ments of the equation of state in the high-temperature 
regime might reveal if 65 changes sign as a function 
of temperature. More generally, we find that the low- 
temperature limit of bn at unitarity is fully determined 
by the low-temperature behavior of 62 and Qi. To ar- 
rive at this result, we derive explicit expressions for A6„ 
and 65^°^ for n < 20, and determine the low-temperature 
behavior of all terms that enter into A6„ and b^f. Us- 
ing the ground state energies at unitarity for trapped 
two-component Fermi gases with up to n = 20 (with 
|ni— 7X2 1 = or 1) [1] , we find that A6„ falls off faster than 
b'i^^ with decreasing T, thus allowing us to obtain ana- 
lytic expressions for the leading order low-temperature 
behavior of &„: b2n exp[i7r(n — l)]exp[— (2n — 
3/2)a)]/(2n) and &2n+i exp(i7rri) exp(— 2na;) for 

n = 1,2,- ••. Thus, the sign of 6„ in the low- 
temperature regime is —,—,+,+,—,—,+, •• • for 
n = 2, 3, 4, 5, 6, 7, 8, 9, 10, • • • . For ri = 2 - 5, we have 
checked that these analytical predictions agree with our 
numerically determined virial coefficients in the low- 
temperature regime. While the sign of 6„ in the low- 
temperature regime may not allow one to draw conclu- 
sions about bl^\ it is interesting, at least from a theoret- 
ical point of view, that the sign and functional form of 
bn in the low-temperaure regime are fully determined by 
Qi and &2- 



V. SUMMARY 

This paper considered the energy spectra of small 
trapped two-component Fermi gases with vanishing and 
finite angular momentum as well as natural and unnatu- 
ral parity. Large portions of the energy spectra of the 
(2,2) and (3,1) systems at unitarity were determined 
as a function of the range of the underlying two-body 
model potential and extrapolated to the zero-range limit. 
The extrapolated zero-range energies are expected to be 
universal, i.e., independent of the underlying Gaussian 
model potential. Portions of the energy spectra of the 
(3, 2) and (4, 1) systems at unitarity were also deter- 
mined. The energies were obtained by solving the rela- 
tive Schrodinger equation using the stochastic variational 
approach. Compact expressions for the relevant matrix 
elements were presented in the appendix. 

The (2, 2) and (3, 1) energies at unitarity were then 
used to determine the low-temperature behavior of the 
fourth-order virial coefficient 64 of the trapped Fermi gas. 
The high-temperature limit of the fourth-order virial co- 
efficient enters into the universal virial equation of state. 
The present study suggests that much larger portions of 
the microscopic energy spectra are needed to predict the 
high-temperature limit of 64. In our view this is unfortu- 
nate. Despite this, we believe that the analysis presented 
illuminates important characteristics relevant to the de- 
termination of the virial coefficients. 

We gratefully acknowledge support by the ARO. KMD 
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Appendix A: Matrix elements 

This appendix summarizes the expressions for the over- 
lap, kinetic energy, trap potential, and interaction poten- 
tial matrix elements for states with natural parity (any 
L) and unnatural parity (L > 0). For notational sim- 
plicity, we omit the subscripts of the matrix A^. and the 
vectors uik and U2k, and consider the matrix elements be- 
tween the unsymmetrized basis functions ip and ip' char- 
acterized by {A,ui,U2) and {A' ,u'i,U2), respectively [see 
Eq. ([5]) of Sec. [n]. The matrix elements have been de- 
rived in the literature [33 - l37| and are summarized here 
for completeness. 

Before providing explicit expressions for the matrix el- 
ements, we introduce a number of auxiliary quantities 
that are utilized in Subsecs. lA II and lA 21 The product 
of ^' and Ip can be conveniently written in terms of the 
matrix B_, 

B^A! + A. (Al) 
We further define the scalars C and pij = 1 or 2), 

and 

Prj = i<fB-%- (A3) 

note that the order of the primed and unprimed vectors 
and Uj matters. We further define the scalars R and 
Sij {hi = 1 or 2), 

R = 3Tr(B"^AAA') (A4) 

and 

= {u,fB-^AMB-^Uj, (A5) 

where the diagonal elements of the matrix A are given 
by the inverse of the masses associated with the Jacobi 
vectors and the off-diagonal elements of A are zero. In 
Eq. (|A4[) . Tr denotes the trace operator. The scalars 

and (p = 1, • • • , n and q = p + 1, • • • , n) have 

a similar structure to R and Sij , 

i?(P9) = 3Tr(B-iQ(f«)) (A6) 

and 

^Ir' - {u[fB-^Q}"'^B-%. (A7) 



The matrix Q^^'^' is defined as 

Q^Pi) ^ ^(Pi) {ujiPi^Y , (A8) 

where a;(P«) is the {n — l)-dimensional vector that re- 
lates the distance vectors r^q to the Jacobi vectors x = 
(pi, • • • ,p„-i), 




(A9) 



Lastly, we define the total mass Mtot, 

Mtot = mp. (AlO) 



1. Natural parity 

For natural parity states, we use h = L and ^2 = in 
Eq. ([5]), which implies that ip' and ip are independent of 
1*2 and U2, respectively. In the following, we assume that 
Ip' and ip are characterizd by the same L and 11 values. 
Under these assumptions the overlap matrix element is 
given by 

{ij'm)^Nr'cpi,, (All) 

where N'l'^^ is a L-dependent constant that enters into 
all matrix elements and thus cancels when calculating 
expectation values. The kinetic energy matrix element 
reads 

(V^'lr-^IV) = iVr*y C(i?pii + 2L5n)pn. (A12) 
The matrix element for the trapping potential reads 

p=l,q>p ^ ^ 

(i?(f'Vii + 2i5ir^)pii.(A13) 

Lastly, the interaction matrix element for the Gaussian 
potential can be written as 

ni n 

-VoY. E {^'\eM-rlq/{2rl)m. (A14) 

p=l q=ni+l 

The expression for the matrix element 
{ip'\exp[—rpg/{2rQ)]\'ilj) reduces to that for the overlap 
matrix element if the matrices -A' and A are replaced by 
A' + Q(f'V(2»-o) and A + Q^p"'' /(2r2), respectively. 
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2. Unnatural parity (L > 0) 



The matrix element for the trapping potential reads 



For unnatural parity states with L > 0, we use li ^ L 
and ?2 = 1 in Eq. ([5]). In the following, we assume that 
tfj' and -0 are characterizd by the same L and 11 values. 
Under these assumptions the overlap matrix element is 
given by 

(ij'm) = iV— Cpii(piip22 " P12P21), (A15) 

where iV]J""^' is a L-dependent constant that enters into 
all matrix elements and thus cancels when calculating 
expectation values. The kinetic energy matrix element 
reads 



{ 

2pii 



p=l,q>p 
c(P'3)l 



2(L-l)5ir (P11P22 - P12P21) 



P11D22 



P22^ir^-pi25r-p21^ir)XAl7) 



[pi) 



{[Rpii + 2(i - l)5'ii] (pii/922 - P12P21) + 

2pil (Pll<5'22 + P22S11 — P12S2I — P21<5'l2)}- 



(A16) 



As in the natural parity case, the expression for the inter- 
action matrix element for the Gaussian potential can be 
related to that of the overlap matrix element by making 
the appropriate substitutions. 
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